Sector-coupling with pypsa#

Note

If you have not yet set up Python on your computer, you can execute this tutorial in your browser via Google Colab. Click on the rocket in the top right corner and launch “Colab”. If that doesn’t work download the .ipynb file and import it in Google Colab.

Then install the following packages by executing the following command in a Jupyter cell at the top of the notebook.

!pip install pypsa pandas numpy matplotlib plotly highspy
import pypsa
import pandas as pd
import numpy as np
import plotly.io as pio
import plotly.offline as py
import matplotlib.pyplot as plt
pd.options.plotting.backend = "plotly"

Previous electricity-only PyPSA model#

To explore sector-coupling options with PyPSA, let’s load the capacity expansion model we built for the electricity system and add sector-coupling technologies and demands on top.

Some of the sector-coupling technologies we are going to add have multiple ouputs (e.g. CHP plants producing heat and power). PyPSA can automatically handle links have more than one input (bus0) and/or output (i.e. bus1, bus2, bus3) with a given efficieny (efficiency, efficiency2, efficiency3).

url = "https://tubcloud.tu-berlin.de/s/Rz4iDEQq8DwBd6R/download/electricity-network.nc"
n = pypsa.Network(url)
INFO:pypsa.io:Retrieving network data from https://tubcloud.tu-berlin.de/s/Rz4iDEQq8DwBd6R/download/electricity-network.nc
INFO:pypsa.io:Imported network electricity-network.nc has buses, carriers, generators, global_constraints, loads, storage_units
n
PyPSA Network
Components:
 - Bus: 1
 - Carrier: 6
 - Generator: 4
 - GlobalConstraint: 1
 - Load: 1
 - StorageUnit: 2
Snapshots: 2190
n.buses.index
Index(['electricity'], dtype='object', name='Bus')
n.generators.index
Index(['OCGT', 'onwind', 'offwind', 'solar'], dtype='object', name='Generator')
n.storage_units.index
Index(['battery storage', 'hydrogen storage underground'], dtype='object', name='StorageUnit')

Hydrogen Production#

The following example shows how to model the components of hydrogen storage separately, i.e. electrolysis, fuel cell and storage.

First, let’s remove the simplified hydrogen storage representation:

n.remove("StorageUnit", "hydrogen storage underground")

Add a separate Bus for the hydrogen energy carrier:

n.add("Bus", "hydrogen", carrier='hydrogen')
Index(['hydrogen'], dtype='object')

Add a Link for the hydrogen electrolysis:

n.add(
    "Link",
    "electrolysis",
    bus0="electricity",
    bus1="hydrogen",
    carrier="electrolysis",
    p_nom_extendable=True,
    efficiency=0.7,
    capital_cost=50e3,  # €/MW/a
)
Index(['electrolysis'], dtype='object')

Add a Link for the fuel cell which reconverts hydrogen to electricity:

n.add(
    "Link",
    "fuel cell",
    bus0="hydrogen",
    bus1="electricity",
    carrier="fuel cell",
    p_nom_extendable=True,
    efficiency=0.5,
    capital_cost=120e3,  # €/MW/a
)
Index(['fuel cell'], dtype='object')

Add a Store for the hydrogen storage:

n.add(
    "Store",
    "hydrogen storage",
    bus="hydrogen",
    carrier="hydrogen storage",
    capital_cost=140,  # €/MWh/a
    e_nom_extendable=True,
    e_cyclic=True,  # cyclic state of charge
)
Index(['hydrogen storage'], dtype='object')

We can also add a hydrogen demand to the hydrogen bus.

In the example below, we add a constant hydrogen demand the size of the electricity demand.

p_set = n.loads_t.p_set["demand"].mean()
p_set
54671.88812785388
n.add("Load", "hydrogen demand", bus="hydrogen", carrier="hydrogen", p_set=p_set)  # MW
Index(['hydrogen demand'], dtype='object')

When we now optimize the model with additional hydrogen demand…

n.optimize()
WARNING:pypsa.consistency:The following links have carriers which are not defined:
Index(['electrolysis', 'fuel cell'], dtype='object', name='Link')
WARNING:pypsa.consistency:The following stores have carriers which are not defined:
Index(['hydrogen storage'], dtype='object', name='Store')
WARNING:pypsa.consistency:The following loads have carriers which are not defined:
Index(['hydrogen demand'], dtype='object', name='Load')
WARNING:pypsa.consistency:The following buses have carriers which are not defined:
Index(['electricity', 'hydrogen'], dtype='object', name='Bus')
WARNING:pypsa.consistency:The following links have carriers which are not defined:
Index(['electrolysis', 'fuel cell'], dtype='object', name='Link')
WARNING:pypsa.consistency:The following stores have carriers which are not defined:
Index(['hydrogen storage'], dtype='object', name='Store')
WARNING:pypsa.consistency:The following loads have carriers which are not defined:
Index(['hydrogen demand'], dtype='object', name='Load')
WARNING:pypsa.consistency:The following buses have carriers which are not defined:
Index(['electricity', 'hydrogen'], dtype='object', name='Bus')
/opt/hostedtoolcache/Python/3.12.7/x64/lib/python3.12/site-packages/linopy/common.py:147: UserWarning:

coords for dimension(s) ['Generator'] is not aligned with the pandas object. Previously, the indexes of the pandas were ignored and overwritten in these cases. Now, the pandas object's coordinates are taken considered for alignment.

INFO:linopy.model: Solve problem using Highs solver
INFO:linopy.io:Writing objective.
Writing constraints.:   0%|          | 0/26 [00:00<?, ?it/s]
Writing constraints.:  38%|███▊      | 10/26 [00:00<00:00, 92.71it/s]
Writing constraints.:  77%|███████▋  | 20/26 [00:00<00:00, 76.21it/s]
Writing constraints.: 100%|██████████| 26/26 [00:00<00:00, 66.08it/s]

Writing continuous variables.:   0%|          | 0/11 [00:00<?, ?it/s]
Writing continuous variables.: 100%|██████████| 11/11 [00:00<00:00, 160.08it/s]
INFO:linopy.io: Writing time: 0.49s
INFO:linopy.solvers:Log file at /tmp/highs.log
Running HiGHS 1.8.0 (git hash: eda5cbe): Copyright (c) 2024 HiGHS under MIT licence terms
Coefficient ranges:
  Matrix [1e-04, 6e+00]
  Cost   [4e-02, 2e+05]
  Bound  [0e+00, 0e+00]
  RHS    [3e+04, 8e+04]
Presolving model
27420 rows, 20857 cols, 73500 nonzeros  0s
25230 rows, 18667 cols, 69120 nonzeros  0s
Presolve : Reductions: rows 25230(-27343); columns 18667(-5431); elements 69120(-45908)
Solving the presolved LP
Using EKK dual simplex solver - serial
  Iteration        Objective     Infeasibilities num(sum)
          0     0.0000000000e+00 Pr: 4380(8.00224e+09) 0s
      17435     7.0371994395e+10 Pr: 8886(5.10211e+11); Du: 0(1.21526e-08) 5s
INFO:linopy.constants: Optimization successful: 
Status: ok
Termination condition: optimal
Solution: 24098 primals, 52573 duals
Objective: 9.45e+10
Solver model: available
Solver message: optimal
      26444     9.4542662147e+10 Pr: 0(0); Du: 0(2.40086e-14) 10s
Solving the original LP from the solution after postsolve
Model   status      : Optimal
Simplex   iterations: 26444
Objective value     :  9.4542662147e+10
HiGHS run time      :          9.96
Writing the solution to /tmp/linopy-solve-dcjh71u4.sol
INFO:pypsa.optimization.optimize:The shadow-prices of the constraints Generator-ext-p-lower, Generator-ext-p-upper, Link-ext-p-lower, Link-ext-p-upper, Store-ext-e-lower, Store-ext-e-upper, StorageUnit-ext-p_dispatch-lower, StorageUnit-ext-p_dispatch-upper, StorageUnit-ext-p_store-lower, StorageUnit-ext-p_store-upper, StorageUnit-ext-state_of_charge-lower, StorageUnit-ext-state_of_charge-upper, StorageUnit-energy_balance, Store-energy_balance were not assigned to the network.
('ok', 'optimal')

… we get an additional price time series for the hydrogen “market”,

n.buses_t.marginal_price.plot()

and can see the individual sizing of the electrolyser, fuel cell and hydrogen storage:

n.statistics.expanded_capacity().div(1e3).round(1)
component    carrier         
Link         electrolysis          183.3
             fuel cell              71.3
Store        hydrogen storage    41979.1
StorageUnit  battery storage        65.1
Generator    offwind               179.2
             onwind                 80.3
             solar                 480.3
dtype: float64

Furthermore, we might want to explore the storage state of charge of the hydrogen storage and the balancing patterns:

n.stores_t.e.plot()

Heat Demand#

For modelling simple heating systems, we create another bus and connect a load with the heat demand time series to it:

n.add("Bus", "heat", carrier='heat')
Index(['heat'], dtype='object')
url = "https://tubcloud.tu-berlin.de/s/mSkHERH8fJCKNXx/download/heat-load-example.csv"
p_set = pd.read_csv(url, index_col=0, parse_dates=True).squeeze()
p_set.head()
snapshot
2015-01-01 00:00:00     61726.043437
2015-01-01 04:00:00    108787.133591
2015-01-01 08:00:00    101508.988082
2015-01-01 12:00:00     90475.260586
2015-01-01 16:00:00     96307.755312
Name: 0, dtype: float64
n.add("Load", "heat demand", carrier="heat", bus="heat", p_set=p_set)
Index(['heat demand'], dtype='object')
n.loads_t.p_set.div(1e3).plot()

What is now missing are a few heat supply options…

Heat pumps#

To model heat pumps, first we have to calculate the coefficient of performance (COP) profile based on the temperature profile of the heat source.

In the example below, we calculate the COP for an air-sourced heat pump with a sink temperature of 55° C and a population-weighted ambient temperature profile for Germany.

The heat pump performance is assumed to be given by the following function:

\[ COP(\Delta T) = 6.81 - 0.121 \Delta T + 0.00063^\Delta T^2 \]

where \(\Delta T = T_{sink} - T_{source}\).

def cop(t_source, t_sink=55):
    delta_t = t_sink - t_source
    return 6.81 - 0.121 * delta_t + 0.000630 * delta_t**2
url = "https://tubcloud.tu-berlin.de/s/S4jRAQMP5Te96jW/download/ninja_weather_country_DE_merra-2_population_weighted.csv"
temp = pd.read_csv(url, skiprows=2, index_col=0, parse_dates=True).loc[
    "2015", "temperature"
][::4]
cop(temp).plot()
plt.scatter(temp, cop(temp))
plt.xlabel("temperature [°C]")
plt.ylabel("COP [-]")
Text(0, 0.5, 'COP [-]')
_images/a06555c6b6757122e19fa62f134d748d0f8d7382acd4797cfe8eff0e3fc90cb5.png

Once we have calculated the heat pump coefficient of performance, we can add the heat pump to the network as a Link. We use the parameter efficiency to incorporate the COP.

n.add(
    "Link",
    "heat pump",
    carrier="heat pump",
    bus0="electricity",
    bus1="heat",
    efficiency=cop(temp),
    p_nom_extendable=True,
    capital_cost=3e5,  # €/MWe/a
)
Index(['heat pump'], dtype='object')

Let’s also add a resistive heater as backup technology:

n.add(
    "Link",
    "resistive heater",
    carrier="resistive heater",
    bus0="electricity",
    bus1="heat",
    efficiency=0.9,
    capital_cost=1e4,  # €/MWe/a
    p_nom_extendable=True,
)
Index(['resistive heater'], dtype='object')

Combined Heat-and-Power (CHP)#

In the following, we are going to add gas-fired combined heat-and-power plants (CHPs). Today, these would use fossil gas, but in the example below we assume imported green methane or methanol with relatively high marginal costs. Since we have no other net emission technology, we can remove the CO\(_2\) limit.

n.remove("GlobalConstraint", "CO2Limit")

Then, we explicitly represent the energy carrier gas:

n.add("Bus", "gas", carrier="gas")
Index(['gas'], dtype='object')

And add a Store of gas, which can be depleted (up to 100 TWh) with fuel costs of 180 €/MWh.

n.add(
    "Store",
    "gas storage",
    carrier="gas storage",
    e_initial=100e6,  # MWh
    e_nom=100e6,  # MWh
    bus="gas",
    marginal_cost=180,  # €/MWh_th
)
Index(['gas storage'], dtype='object')

When we do this, we have to model the OCGT power plant as link which converts gas to electricity, not as generator.

n.remove("Generator", "OCGT")
n.add(
    "Link",
    "OCGT",
    bus0="gas",
    bus1="electricity",
    carrier="OCGT",
    p_nom_extendable=True,
    capital_cost=20000,  # €/MW/a
    efficiency=0.4,
)
Index(['OCGT'], dtype='object')

Next, we are going to add a combined heat-and-power (CHP) plant with fixed heat-power ratio (i.e. backpressure operation).

Note

If you want to model flexible heat-power ratios, have a look at this example: https://pypsa.readthedocs.io/en/latest/examples/power-to-gas-boiler-chp.html

n.add(
    "Link",
    "CHP",
    bus0="gas",
    bus1="electricity",
    bus2="heat",
    carrier="CHP",
    p_nom_extendable=True,
    capital_cost=40000,
    efficiency=0.4,
    efficiency2=0.4,
)
Index(['CHP'], dtype='object')

Now, let’s optimize the current status of model:

n.optimize()
WARNING:pypsa.consistency:The following links have carriers which are not defined:
Index(['electrolysis', 'fuel cell', 'heat pump', 'resistive heater', 'CHP'], dtype='object', name='Link')
WARNING:pypsa.consistency:Encountered nan's in static data for columns ['efficiency2'] of component 'Link'.
WARNING:pypsa.consistency:The following stores have carriers which are not defined:
Index(['hydrogen storage', 'gas storage'], dtype='object', name='Store')
WARNING:pypsa.consistency:The following loads have carriers which are not defined:
Index(['hydrogen demand', 'heat demand'], dtype='object', name='Load')
WARNING:pypsa.consistency:The following buses have carriers which are not defined:
Index(['electricity', 'hydrogen', 'heat', 'gas'], dtype='object', name='Bus')
WARNING:pypsa.consistency:The following links have carriers which are not defined:
Index(['electrolysis', 'fuel cell', 'heat pump', 'resistive heater', 'CHP'], dtype='object', name='Link')
WARNING:pypsa.consistency:Encountered nan's in static data for columns ['efficiency2'] of component 'Link'.
WARNING:pypsa.consistency:The following stores have carriers which are not defined:
Index(['hydrogen storage', 'gas storage'], dtype='object', name='Store')
WARNING:pypsa.consistency:The following loads have carriers which are not defined:
Index(['hydrogen demand', 'heat demand'], dtype='object', name='Load')
WARNING:pypsa.consistency:The following buses have carriers which are not defined:
Index(['electricity', 'hydrogen', 'heat', 'gas'], dtype='object', name='Bus')
/opt/hostedtoolcache/Python/3.12.7/x64/lib/python3.12/site-packages/linopy/common.py:147: UserWarning:

coords for dimension(s) ['Generator'] is not aligned with the pandas object. Previously, the indexes of the pandas were ignored and overwritten in these cases. Now, the pandas object's coordinates are taken considered for alignment.

INFO:linopy.model: Solve problem using Highs solver
INFO:linopy.io:Writing objective.
Writing constraints.:   0%|          | 0/27 [00:00<?, ?it/s]
Writing constraints.:  41%|████      | 11/27 [00:00<00:00, 100.28it/s]
Writing constraints.:  81%|████████▏ | 22/27 [00:00<00:00, 79.38it/s] 
Writing constraints.: 100%|██████████| 27/27 [00:00<00:00, 59.41it/s]

Writing continuous variables.:   0%|          | 0/11 [00:00<?, ?it/s]
Writing continuous variables.: 100%|██████████| 11/11 [00:00<00:00, 150.88it/s]
INFO:linopy.io: Writing time: 0.55s
INFO:linopy.solvers:Log file at /tmp/highs.log
Running HiGHS 1.8.0 (git hash: eda5cbe): Copyright (c) 2024 HiGHS under MIT licence terms
Coefficient ranges:
  Matrix [1e-04, 6e+00]
  Cost   [4e-02, 3e+05]
  Bound  [0e+00, 0e+00]
  RHS    [2e+04, 1e+08]
Presolving model
42749 rows, 34000 cols, 119487 nonzeros  0s
38370 rows, 29621 cols, 110729 nonzeros  0s
Presolve : Reductions: rows 38370(-38294); columns 29621(-5430); elements 110729(-50291)
Solving the presolved LP
Using EKK dual simplex solver - serial
  Iteration        Objective     Infeasibilities num(sum)
          0     0.0000000000e+00 Ph1: 0(0) 0s
      15570     2.2962917512e+10 Pr: 7061(9.5299e+11); Du: 0(1.51356e-06) 5s
      21609     3.7705371248e+10 Pr: 8310(1.11402e+13); Du: 0(2.50085e-07) 10s
      26715     9.9340691224e+10 Pr: 12410(2.58287e+13); Du: 0(7.68579e-08) 16s
      31162     1.3762022038e+11 Pr: 15176(1.90995e+12) 22s
      33878     1.4301664789e+11 Pr: 0(0); Du: 0(3.73035e-14) 25s
Solving the original LP from the solution after postsolve
Model   status      : Optimal
Simplex   iterations: 33878
Objective value     :  1.4301664789e+11
HiGHS run time      :         24.83
INFO:linopy.constants: Optimization successful: 
Status: ok
Termination condition: optimal
Solution: 35051 primals, 76664 duals
Objective: 1.43e+11
Solver model: available
Solver message: optimal
INFO:pypsa.optimization.optimize:The shadow-prices of the constraints Generator-ext-p-lower, Generator-ext-p-upper, Link-ext-p-lower, Link-ext-p-upper, Store-fix-e-lower, Store-fix-e-upper, Store-ext-e-lower, Store-ext-e-upper, StorageUnit-ext-p_dispatch-lower, StorageUnit-ext-p_dispatch-upper, StorageUnit-ext-p_store-lower, StorageUnit-ext-p_store-upper, StorageUnit-ext-state_of_charge-lower, StorageUnit-ext-state_of_charge-upper, StorageUnit-energy_balance, Store-energy_balance were not assigned to the network.
Writing the solution to /tmp/linopy-solve-wn7t42yh.sol
('ok', 'optimal')

The objective cost in bn€/a:

n.objective / 1e9
143.01664789325199

The heat energy balance (positive is supply, negative is consumption):

n.statistics.energy_balance(bus_carrier='heat').div(1e6).round(1)
component  carrier           bus_carrier
Link       heat pump         heat           519.1
           resistive heater  heat            34.6
           CHP               heat            40.0
Load       heat              heat          -593.7
dtype: float64

The heat energy balance as a time series:

n.statistics.energy_balance(aggregate_time=False, bus_carrier='heat').div(1e3).groupby("carrier").sum().T.plot()

The price time series:

n.buses_t.marginal_price.plot()

Long-duration heat storage#

One technology of particular interest in district heating systems with large shares of renewables is long-duration thermal energy storage.

In the following, we are going to introduce a heat storage with investment cost of approximately 3 €/kWh. The energy is not perfectly stored in water tanks. There are standing losses. The decay of thermal energy in the heat storage is modelled through the function \(1 - e^{-\frac{1}{24\tau}}\), where \(\tau\) is assumed to be 180 days. We want to see how that influences the optimal design decisions in the heating sector.

n.add(
    "Store",
    "heat storage",
    bus="heat",
    carrier="heat storage",
    capital_cost=300, # roughly annuity of 3 €/kWh
    standing_loss=1 - np.exp(-1 / 24 / 180),
    e_nom_extendable=True,
)
Index(['heat storage'], dtype='object')
n.optimize()
WARNING:pypsa.consistency:The following links have carriers which are not defined:
Index(['electrolysis', 'fuel cell', 'heat pump', 'resistive heater', 'CHP'], dtype='object', name='Link')
WARNING:pypsa.consistency:Encountered nan's in static data for columns ['efficiency2'] of component 'Link'.
WARNING:pypsa.consistency:The following stores have carriers which are not defined:
Index(['hydrogen storage', 'gas storage', 'heat storage'], dtype='object', name='Store')
WARNING:pypsa.consistency:The following loads have carriers which are not defined:
Index(['hydrogen demand', 'heat demand'], dtype='object', name='Load')
WARNING:pypsa.consistency:The following buses have carriers which are not defined:
Index(['electricity', 'hydrogen', 'heat', 'gas'], dtype='object', name='Bus')
WARNING:pypsa.consistency:The following links have carriers which are not defined:
Index(['electrolysis', 'fuel cell', 'heat pump', 'resistive heater', 'CHP'], dtype='object', name='Link')
WARNING:pypsa.consistency:Encountered nan's in static data for columns ['efficiency2'] of component 'Link'.
WARNING:pypsa.consistency:The following stores have carriers which are not defined:
Index(['hydrogen storage', 'gas storage', 'heat storage'], dtype='object', name='Store')
WARNING:pypsa.consistency:The following loads have carriers which are not defined:
Index(['hydrogen demand', 'heat demand'], dtype='object', name='Load')
WARNING:pypsa.consistency:The following buses have carriers which are not defined:
Index(['electricity', 'hydrogen', 'heat', 'gas'], dtype='object', name='Bus')
/opt/hostedtoolcache/Python/3.12.7/x64/lib/python3.12/site-packages/linopy/common.py:147: UserWarning:

coords for dimension(s) ['Generator'] is not aligned with the pandas object. Previously, the indexes of the pandas were ignored and overwritten in these cases. Now, the pandas object's coordinates are taken considered for alignment.

INFO:linopy.model: Solve problem using Highs solver
INFO:linopy.io:Writing objective.
Writing constraints.:   0%|          | 0/27 [00:00<?, ?it/s]
Writing constraints.:  41%|████      | 11/27 [00:00<00:00, 101.35it/s]
Writing constraints.:  81%|████████▏ | 22/27 [00:00<00:00, 73.57it/s] 
Writing constraints.: 100%|██████████| 27/27 [00:00<00:00, 59.19it/s]

Writing continuous variables.:   0%|          | 0/11 [00:00<?, ?it/s]
Writing continuous variables.: 100%|██████████| 11/11 [00:00<00:00, 137.22it/s]
INFO:linopy.io: Writing time: 0.56s
INFO:linopy.solvers:Log file at /tmp/highs.log
Running HiGHS 1.8.0 (git hash: eda5cbe): Copyright (c) 2024 HiGHS under MIT licence terms
Coefficient ranges:
  Matrix [1e-04, 6e+00]
  Cost   [4e-02, 3e+05]
  Bound  [0e+00, 0e+00]
  RHS    [2e+04, 1e+08]
Presolving model
47128 rows, 38380 cols, 132624 nonzeros  0s
40560 rows, 31812 cols, 119488 nonzeros  0s
40560 rows, 31812 cols, 119488 nonzeros  0s
Presolve : Reductions: rows 40560(-42675); columns 31812(-7620); elements 119488(-56862)
Solving the presolved LP
Using EKK dual simplex solver - serial
  Iteration        Objective     Infeasibilities num(sum)
          0     0.0000000000e+00 Ph1: 0(0) 0s
      16443     7.6508187661e+09 Pr: 8788(1.48947e+14); Du: 0(9.70147e-07) 5s
      22711     4.3725291538e+10 Pr: 6314(1.40186e+13); Du: 0(1.2653e-06) 11s
      27869     4.4056163819e+10 Pr: 7072(4.98632e+14); Du: 0(6.05973e-07) 16s
      31942     6.7369736899e+10 Pr: 10332(3.3962e+15); Du: 0(2.37077e-07) 21s
      35534     9.5675850589e+10 Pr: 18650(1.00973e+14) 27s

The objective cost (in bn€/a) was reduced by 7%!

n.objective / 1e9
133.41557805470643

The heat energy balance shows the additional losses of the heat storage and the added supply:

n.statistics.energy_balance(bus_carrier='heat').div(1e6).round(1)
component  carrier           bus_carrier
Load       heat              heat          -593.7
Store      heat storage      heat           -11.7
Link       heat pump         heat           499.2
           resistive heater  heat            82.3
           CHP               heat            23.8
dtype: float64

The heat energy balance as a time series:

n.statistics.energy_balance(aggregate_time=False, bus_carrier='heat').div(1e3).groupby("carrier").sum().T.plot()

The heat price time series is much less impacted by the variability of the electricity price:

n.buses_t.marginal_price.plot()

The different storage state of charge time series:

n.stores_t.e.plot()

Electric Vehicles#

To model electric vehicles, we first create another bus for the electric vehicles.

n.add("Bus", "EV", carrier="EV")
Index(['EV'], dtype='object')

Then, we can attach the electricity consumption of electric vehicles to this bus:

url = "https://tubcloud.tu-berlin.de/s/9r5bMSbzzQiqG7H/download/electric-vehicle-profile-example.csv"
p_set = pd.read_csv(url, index_col=0, parse_dates=True).squeeze()
p_set.loc["2015-01"].div(1e3).plot()
n.add("Load", "EV demand", bus="EV", carrier="EV demand", p_set=p_set)
Index(['EV demand'], dtype='object')

The electric vehicles can only be charged when they are plugged-in. Below we load an availability profile telling us what share of electric vehicles is plugged-in at home – we only assume home charging in this example.

url = "https://tubcloud.tu-berlin.de/s/E3PBWPfYaWwCq7a/download/electric-vehicle-availability-example.csv"
availability_profile = pd.read_csv(url, index_col=0, parse_dates=True).squeeze()
availability_profile.loc["2015-01"].plot()

Then, we can add a link for the electric vehicle charger using assumption about the number of EVs and their charging rates.

number_cars = 40e6  #  number of EV cars
bev_charger_rate = 0.011  # 3-phase EV charger with 11 kW
p_nom = number_cars * bev_charger_rate
n.add(
    "Link",
    "EV charger",
    bus0="electricity",
    bus1="EV",
    p_nom=p_nom,
    carrier="EV charger",
    p_max_pu=availability_profile,
    efficiency=0.9,
)
Index(['EV charger'], dtype='object')

We can also allow vehicle-to-grid operation (i.e. electric vehicles inject power into the grid):

n.add(
    "Link",
    "V2G",
    bus0="EV",
    bus1="electricity",
    p_nom=p_nom,
    carrier="V2G",
    p_max_pu=availability_profile,
    efficiency=0.9,
)
Index(['V2G'], dtype='object')

The demand-side management potential we model as a store. This is not unlike a battery storage, but we impose additional constraints on when the store needs to be charged to a certain level (e.g. 75% full every morning).

bev_energy = 0.05  # average battery size of EV in MWh
bev_dsm_participants = 0.5  # share of cars that do smart charging

e_nom = number_cars * bev_energy * bev_dsm_participants
url = "https://tubcloud.tu-berlin.de/s/K62yACBRTrxLTia/download/dsm-profile-example.csv"
dsm_profile = pd.read_csv(url, index_col=0, parse_dates=True).squeeze().shift(2).fillna(0)[::4]
dsm_profile.loc["2015-01-01":"2015-01-07"].plot()
n.add(
    "Store",
    "EV DSM",
    bus="EV",
    carrier="EV battery",
    e_cyclic=True,  # state of charge at beginning = state of charge at the end
    e_nom=e_nom,
    e_min_pu=dsm_profile,
)
Index(['EV DSM'], dtype='object')

Then, we can solve the fully sector-coupled model altogether including electricity, passenger transport, hydrogen and heating.

n.optimize(solver_name="highs")
Hide code cell output
WARNING:pypsa.consistency:The following stores have carriers which are not defined:
Index(['hydrogen storage', 'gas storage', 'heat storage', 'EV DSM'], dtype='object', name='Store')
WARNING:pypsa.consistency:The following buses have carriers which are not defined:
Index(['electricity', 'hydrogen', 'heat', 'gas', 'EV'], dtype='object', name='Bus')
WARNING:pypsa.consistency:The following loads have carriers which are not defined:
Index(['hydrogen demand', 'heat demand', 'EV demand'], dtype='object', name='Load')
WARNING:pypsa.consistency:The following links have carriers which are not defined:
Index(['electrolysis', 'fuel cell', 'heat pump', 'resistive heater', 'CHP',
       'EV charger', 'V2G'],
      dtype='object', name='Link')
WARNING:pypsa.consistency:The following stores have carriers which are not defined:
Index(['hydrogen storage', 'gas storage', 'heat storage', 'EV DSM'], dtype='object', name='Store')
WARNING:pypsa.consistency:The following buses have carriers which are not defined:
Index(['electricity', 'hydrogen', 'heat', 'gas', 'EV'], dtype='object', name='Bus')
WARNING:pypsa.consistency:The following loads have carriers which are not defined:
Index(['hydrogen demand', 'heat demand', 'EV demand'], dtype='object', name='Load')
WARNING:pypsa.consistency:The following links have carriers which are not defined:
Index(['electrolysis', 'fuel cell', 'heat pump', 'resistive heater', 'CHP',
       'EV charger', 'V2G'],
      dtype='object', name='Link')
/home/fneum/miniconda3/envs/nordnet/lib/python3.12/site-packages/linopy/common.py:147: UserWarning:

coords for dimension(s) ['Generator'] is not aligned with the pandas object. Previously, the indexes of the pandas were ignored and overwritten in these cases. Now, the pandas object's coordinates are taken considered for alignment.

INFO:linopy.model: Solve problem using Highs solver
INFO:linopy.io:Writing objective.
Writing constraints.: 100%|██████████| 29/29 [00:00<00:00, 84.81it/s] 
Writing continuous variables.: 100%|██████████| 11/11 [00:00<00:00, 226.05it/s]
INFO:linopy.io: Writing time: 0.41s
INFO:linopy.solvers:Log file at /tmp/highs.log
Running HiGHS 1.7.2 (git hash: 184e327): Copyright (c) 2024 HiGHS under MIT licence terms
Coefficient ranges:
  Matrix [1e-04, 6e+00]
  Cost   [4e-02, 3e+05]
  Bound  [0e+00, 0e+00]
  RHS    [2e+03, 1e+08]
Presolving model
51508 rows, 47140 cols, 150144 nonzeros  0s
42750 rows, 38382 cols, 132628 nonzeros  0s
42750 rows, 38382 cols, 132628 nonzeros  0s
Presolve : Reductions: rows 42750(-58005); columns 38382(-9810); elements 132628(-74382)
Solving the presolved LP
Using EKK dual simplex solver - serial
  Iteration        Objective     Infeasibilities num(sum)
          0     0.0000000000e+00 Ph1: 0(0) 0s
      20489     7.0209292127e+09 Pr: 10189(7.06301e+14); Du: 0(9.85806e-07) 5s
      26917     3.1666331621e+10 Pr: 6630(2.35325e+13); Du: 0(6.23236e-07) 11s
      32263     3.2859044901e+10 Pr: 7443(6.50205e+13); Du: 0(2.8441e-07) 16s
      36620     5.3585544208e+10 Pr: 10737(1.8223e+13); Du: 0(4.035e-07) 21s
      40843     8.8898698525e+10 Pr: 15887(1.50189e+15); Du: 0(7.42276e-07) 26s
      45777     1.0089919205e+11 Pr: 9171(1.94314e+16); Du: 0(5.04967e-07) 32s
      49643     1.1413934646e+11 Pr: 14967(2.15477e+13); Du: 0(5.47685e-07) 37s
      53329     1.2988564883e+11 Pr: 10975(1.23989e+13); Du: 0(6.79277e-07) 43s
      56786     1.3519466178e+11 Pr: 4851(3.46768e+11); Du: 0(5.56497e-07) 48s
      59796     1.3651637401e+11 Pr: 0(0); Du: 0(2.27374e-13) 52s
      59796     1.3651637401e+11 Pr: 0(0); Du: 0(2.27374e-13) 52s
Solving the original LP from the solution after postsolve
Model   status      : Optimal
Simplex   iterations: 59796
Objective value     :  1.3651637401e+11
HiGHS run time      :         52.41
Writing the solution to /tmp/linopy-solve-tk463w7g.sol
INFO:linopy.constants: Optimization successful: 
Status: ok
Termination condition: optimal
Solution: 48192 primals, 100755 duals
Objective: 1.37e+11
Solver model: available
Solver message: optimal

INFO:pypsa.optimization.optimize:The shadow-prices of the constraints Generator-ext-p-lower, Generator-ext-p-upper, Link-fix-p-lower, Link-fix-p-upper, Link-ext-p-lower, Link-ext-p-upper, Store-fix-e-lower, Store-fix-e-upper, Store-ext-e-lower, Store-ext-e-upper, StorageUnit-ext-p_dispatch-lower, StorageUnit-ext-p_dispatch-upper, StorageUnit-ext-p_store-lower, StorageUnit-ext-p_store-upper, StorageUnit-ext-state_of_charge-lower, StorageUnit-ext-state_of_charge-upper, StorageUnit-energy_balance, Store-energy_balance were not assigned to the network.
('ok', 'optimal')

A very general overview is given with a call to n.statistics():

n.statistics()
Optimal Capacity Installed Capacity Supply Withdrawal Energy Balance Transmission Capacity Factor Curtailment Capital Expenditure Operational Expenditure Revenue Market Value
Generator offwind 2.455913e+05 0.0 7.805015e+08 0.000000e+00 7.805015e+08 0.0 0.362791 0.000000e+00 4.286949e+10 1.654663e+07 4.288603e+10 54.946768
onwind 1.575662e+05 0.0 2.659010e+08 0.000000e+00 2.659010e+08 0.0 0.192643 1.808124e+07 1.601567e+10 3.798662e+08 1.639554e+10 61.660318
solar 5.033996e+05 0.0 5.476326e+08 0.000000e+00 5.476326e+08 0.0 0.124186 0.000000e+00 2.584798e+10 5.804905e+06 2.585378e+10 47.210086
Link CHP 1.123002e+05 0.0 5.083544e+07 6.354431e+07 -1.270886e+07 0.0 0.064594 0.000000e+00 4.492008e+09 0.000000e+00 4.492008e+09 88.363693
EV charger 4.400000e+05 440000.0 2.000217e+08 2.222463e+08 -2.222463e+07 0.0 0.057660 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 NaN
OCGT 6.770876e+03 0.0 3.536607e+05 8.841518e+05 -5.304911e+05 0.0 0.014907 0.000000e+00 1.354175e+08 0.000000e+00 1.354175e+08 382.902358
V2G 4.400000e+05 440000.0 6.816967e+07 7.574408e+07 -7.574408e+06 0.0 0.019651 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 NaN
electrolysis 1.579354e+05 0.0 4.789257e+08 6.841796e+08 -2.052539e+08 0.0 0.494523 0.000000e+00 7.896770e+09 0.000000e+00 7.896770e+09 16.488505
heat pump 4.444094e+04 0.0 4.804672e+08 1.924107e+08 2.880565e+08 0.0 0.494244 0.000000e+00 1.333228e+10 0.000000e+00 1.333228e+10 27.748578
resistive heater 9.724606e+04 0.0 9.919236e+07 1.102137e+08 -1.102137e+07 0.0 0.129378 0.000000e+00 9.724606e+08 0.000000e+00 9.724606e+08 9.803785
Load - 0.000000e+00 0.0 0.000000e+00 4.789257e+08 -4.789257e+08 0.0 NaN 0.000000e+00 0.000000e+00 0.000000e+00 -4.218900e+10 NaN
EV demand 0.000000e+00 0.0 0.000000e+00 1.242776e+08 -1.242776e+08 0.0 NaN 0.000000e+00 0.000000e+00 0.000000e+00 -1.104624e+10 NaN
heat 0.000000e+00 0.0 0.000000e+00 5.936692e+08 -5.936692e+08 0.0 NaN 0.000000e+00 0.000000e+00 0.000000e+00 -4.801251e+10 NaN
hydrogen 0.000000e+00 0.0 0.000000e+00 4.789257e+08 -4.789257e+08 0.0 NaN 0.000000e+00 0.000000e+00 0.000000e+00 -3.838917e+10 NaN
Store EV battery 1.000000e+06 1000000.0 1.290614e+08 1.290614e+08 0.000000e+00 0.0 0.726787 0.000000e+00 0.000000e+00 0.000000e+00 3.120540e+09 24.178714
gas storage 1.000000e+08 100000000.0 6.442846e+07 0.000000e+00 6.442846e+07 0.0 0.499883 0.000000e+00 0.000000e+00 1.159712e+10 1.159712e+10 180.000000
heat storage 1.270014e+07 0.0 1.953699e+08 2.067779e+08 -1.140802e+07 0.0 0.443182 0.000000e+00 3.810042e+09 0.000000e+00 3.810042e+09 19.501679
hydrogen storage 6.532086e+07 0.0 2.009489e+08 2.009489e+08 0.000000e+00 0.0 0.565556 0.000000e+00 9.144920e+09 0.000000e+00 9.144920e+09 45.508680

And it can tell you statistics about the capital expenditures:

n.statistics.capex().groupby("carrier").sum().div(1e9).sort_values().dropna().plot.bar()

Exercises#

Explore how the model reacts to changing assumptions and available technologies. Here are a few inspirations, but choose in any order according to your interests:

  • By how much did the optimised capacities react to

  • Assume underground hydrogen storage is not geographically available. Increase the cost of hydrogen storage by factor 10. How does the model react?

  • Add a ground-sourced heat pump with a constant COP function of 3.5 but double the investment costs. Would this technology get built? How low would the costs need to be?

  • Vary the charging rates and availability profiles of electric vehicles, and whether they participate in demand-side management (V2G, smart charging) or not. How does it affect the system design?

  • Limit green gas imports to 10 TWh. What does the model do in periods with persistent low wind and solar feed-in but high heating demand?

You can also combine these aspects with investigations suggested for the electricity-only model.